;; A non-robust implementation of the Optimal Bitonic Tour
;; Copyright (C) 2022  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.compgeom.optimal-bitonic-tour)

(defn- dist [[p1 p2] [q1 q2]]
  (let [dx (- q1 p1)
        dy (- q2 p2)]
    (Math/sqrt (+ (* dx dx) (* dy dy)))))

(defn- rule3-reducer [[k-min l-min :as val-min] [k l :as val]]
  (if (< l l-min) val val-min))

(defn- rule3 [front pⱼ L]
  (let [L3 (map #(vector %1 (+ %2 (dist %3 pⱼ))) (range) L front)]
    (reduce rule3-reducer L3)))

;; Reference: https://www.cs.huji.ac.il/course/2004/algo/Solutions/bitonic.pdf
;; (BUT it seems to have been taken down as of 2022 October 16).
;;
;; (L i j) is the length of the shortest normal bitonic path with endpoints
;; pᵢ and pⱼ, 0 ≤ i < j < n:
;; 
;; (L i j) = (dist pᵢ pⱼ)                        if i = 0 and j = 1        (1)
;;         = (+ (L i (dec j)) (dist pⱼ-₁ pⱼ))    if i < j − 1              (2)
;;         = (min₀≤ₖ<ᵢ (+ (L k i) (dist pₖ pⱼ))) if j > 1 and i = j − 1    (3)
;;
;; The algorithm builds L column by column from left to right but only the
;; previous column is needed to compute the next column so we only keep track
;; of the most recent column.
;;
;; N keeps track of which index k resulted in the value of (3).
;; N is indexed by 0 < j < n since i is implied.

(defn- optimal-bitonic-tour-sorted
  "`coll` is assumed to be sorted."
  [coll]
  (if (nil? (second coll))
    [coll nil]
    ; compute the costs for all bitonic tours
    (loop [front (vec (take 2 coll))
           back (seq (drop 2 coll))
           L (vector (dist (first coll) (second coll))) ; (1)
           N [nil 0]]                                   ; (N 0) is unused
      (if back
        (let [pⱼ-₁ (peek front)
              pⱼ (first back)
              d (dist pⱼ-₁ pⱼ)
              L2 (mapv #(+ % d) L)      ; (2)
              [k l] (rule3 front pⱼ L)] ; (3)
          (recur (conj front pⱼ) (next back) (conj L2 l) (conj N k)))
        ; construct the two monotone chains
        (loop [j (dec (count front))
               i (dec j)
               Sj (vector (front j))
               Si (vector (front i))]
          (if (== i 0)
            [(rseq Si) Sj]
            (let [k (N j)]
              (recur (inc k)
                     k
                     (into Si (map front (range (dec i) k -1)))
                     (conj Sj (front k))))))))))

(defn optimal-bitonic-tour
  "Return the shortest bitonic tour of the points in `coll`.
  The return value is a sequence of two monotonic chains.
  The first monotonic chain has non-decreasing x values;
  the second monotonic chain has non-increasing x values."
  [coll]
  (-> coll sort optimal-bitonic-tour-sorted))
